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1. Introduction 

The Space-Time Conservation Element and Solution Element Method, or CE/SE 
Method for short [1-3], is a recently developed numerical method for conservation 
laws. Despite its second order accuracy in space and time, it possesses low dis- 
persion errors and low dissipation. The method is robust enough to cover a wide 
range of compressible flows: from weak linear acoustic waves to strong discontin- 
uous waves (shocks). An outstanding feature of the CE/SE scheme is its truly 
multi-dimensional, simple but effective non-reflecting boundary condition (NRBC), 
which is particularly valuable for CAA (computational aeroacoustics). 

In nature, the method may be categorized as a finite volume method, where 
the conservation element (CE) is equivalent to a finite control volume (or cell) and 
the solution element ( SE ) can be understood as the cell interface. However, due 
to its careful treatment of the surface fluxes and geometry, it is different from the 
existing schemes. 

Currently, the CE/SE scheme has been developed to a matured stage that a 
3-D unstructured CE/SE Navier-Stokes solver is already available. However, in the 
present review paper, as a general introduction to the CE/SE method, only the 2-D 
unstructured Euler CE/SE solver is chosen and sketched in §2. Then applications of 
the 2-D and 3-D CE/SE schemes to linear, and in particular, nonlinear aeroacoustics 
are depicted in § 3, § 4, and § 5 to demonstrate its robustness and capability. 

2. The 2-D CE/SE Euler Scheme 

In this Section, only a brief sketch of the CE/SE scheme for the 2-D Euler 
equations for gas dynamics is given. Detailed description of the method may be 
found in the original papers [1- 3]. 

2.1. Conservation Form of the 2-D Unsteady Euler Equations. Con- 
sider a dimensionless conservation form of the unsteady 2-D Euler equations of a 
perfect gas. Let p, p, and 7 be the density, streamwise and transverse velocity 
components, static pressure, and constant specific heat ratio, respectively. The 2-D 
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hexagon cylinder ADBECFA - A’D’B’E’C’F’A’ and 
its 3 CEs, OADB-O’A’D’B’, OBEC-O’B’E’C’, 
OCFA - O’C’F’A’ 


Figure 1. CE/SE unstructured grid in space-time E 3 space. 

Euler equations then can be written in the following vector form: 

(1) U t + F x + G y = 0, 

where #, y, and t are the streamwise and transversal coordinates and time, re- 
spectively. The conservative flow variable vector U and the flux vectors in the 
streamwise and radial directions, F and G, are given by: 
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with the conservative flow variables 

Ui=p, U 2 = pu, u 3 =pv , U 4 =p/('Y~ 1 ) + p(v? +v 2 )/2. 

where the fluxes are: 

Fi = U 2 , F 2 = ( 7 - 1)U A + [(3 - 7 )Ui - (7 - 1 Wl] /2U U 
F 3 = u 2 u 3 /lh, F 4 = 'yU 2 Ui/U 1 - (7 - 1 )U 2 [U* + U$] /2U?, 

G\ = U 3 , G 2 = U 2 U 3 /lh, G 3 = (j — 1)U 4 + [(3 — 7 )t/f — (7 — 1 )^ 1 ] /2C/i, 
G 4 = jUsUi/lh - (7 - 1 )U 3 [C7| + U$] /2Ul 
By considering (x,y,t) as coordinates of a three-dimensional Euclidean space, E 3 , 
and using Gauss’ divergence theorem, it follows that Eq. ( 1 ) is equivalent to the 
following integral conservation law: 

(2) (f H m • dS = 0, m = 1, 2, 3, 4, 

Js(v ) 

where S(V) denotes the surface around a volume V in E 3 and H m = (F m , G m , U m ). 

In general, the CE/SE method systematically the integral equations (2) and 
naturally captures shocks and other discontinuities in the flow. Both conservative 
variables U and their spatial derivatives TJ X , \J y are computed simultaneously as 
unknowns. 
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2.2. Unstructured Grid for CE/SE. The CE/SE scheme is constructed 
to take advantage of an unstructured triangular grid. The unstructured geometry 
used with the CE/SE scheme is illustrated in Fig. 1. Here, A ABC is a typical 
triangle cell and D,E,F are the triangle centers of the neighboring cells. The 
flow variables U, U®, and at the previous time step are stored at the triangle 
cell centers. Three quadrilateral cylinders (conservation elements) are formed by 
the edges that connect the vertices and the center of the triangle and its three 
neighbors. In the space-time E 3 space, Eq. (2) is applied to the hexagon cylinder 
ADBECF — A'D'B'E'C'F 1 of volume V that consists of these 3 quadrilateral 
cylinder CE$ (Fig 1). 

In the CE/SE scheme, the above flux conservation relation, Eq. (2), in space- 
time is the only mechanism that transfers information between node points. A 
conservation element CE (here, quadrilateral cylinders) is the finite volume to which 
Eq. (2) is applied. Discontinuities are allowed to occur in a conservation element. 
A solution element SE associated with a grid node ( e.g ., D, E, or F in Fig. 1) is 
here a set of interface planes in E 3 that passes through this node (e.g. DAA'D ', 
DBB'D ', EBB'E ', ECC'E ', etc.). Each surface S(CE) is made up of segments 
belonging to two neighboring CE’s. Within a given solution element SE(j, n), 
where j, n are the node index and time step, respectively, the flow variables are not 
only considered continuous but are also approximated by the linear Taylor series 
expansions. The surface flux can thus be calculated accurately and easily by first 
evaluating the flux vectors at the geometrical center of the surface through the 
above Taylor series expansions. 

At time level n, the solution variables U, U x , and U y are given at the three 
nodes D, E, F in Fig. 1 and U, U x and U y at O' at the new time level n + 1 are 
to be evaluated. In principle, each of the 3 CEs provides 4 scalar equations when 
Eq. (2) is applied to the element. Hence, the 12 scalar equations needed for the 
12 scalar unknowns at O' are available. All the unknowns are computed based 
on these relations. No extrapolations (interpolations) across a stencil of cells are 
needed. 


2.3. Non- Reflecting Boundary Conditions (NRBCs) for CE/SE. In 

the CE/SE scheme, NRBCs are constructed so as to allow fluxes remain continuous 
across the boundary surfaces. For example, at the downstream boundary, where 
there are substantial gradients in y direction, the NRBC requires that 

(u x )" = 0, u" = u y 1/2 (Uy)” = (U y )”- 1/2 , 

where j' is the index of an interior node closest to the boundary ghost node j and 
U™ and (U y)™ are now defined by simple extrapolation from the interior. This 
NRBC is valid for either supersonic or subsonic flows. 

In the following sections, the 2-D and 3-D CE/SE Euler (Navier-Stokes) schemes 
are tested to demonstrate its capability and robustness for aeroacoustics computa- 
tions. Several selected problems in linear and nonlinear aeroacoustics computations 
are presented. The numerical results, which cover a wide spectrum of waves, from 
linear and nonlinear acoustic waves to discontinuous waves (shocks), are then com- 
pared to available exact solutions or experimental findings. 
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3. 2-D Linear Aeroacoustics Test Problems 

Three typical linear aeroacoustics examples are considered in this section. More 
details can be found in [5, 6]. 

3.1. Acoustic Pulse, Entropy Wave, and Vorticity Wave Propagation. 

This problem [5] is one of the benchmark problems of the first CAA Workshop 
(Category 3, Problem 1) [4]. The computational domain in the x-y plane is a 
square with —100 < x < 100, and —100 < y < 100. A uniform 201 x 201 grid 
is used with Ax = Ay = 1 . Initially, a Gaussian acoustic pulse is located at the 
center of the domain (0,0) and a weaker entropy /vorticity disturbance is located off 
center (67,0), with a mean flow of Mach number M = 0.5 in the x direction. The 
corresponding non-dimensional pressure, density, and streamwise and transverse 
velocity components are given by 

p = - + Se-^+y 2 ') , p = 1 + <$[ e -“ 1 <* a+ «' 2 ) + 0.1 Se- a ^ x - 67 ^ 2+ y% 

7 

u = M + 0mSye- a2 ^ x - 67 ^ +v2 \ v = -0.04£(> - 67)e- a ^ x - 6 ^ 2+y2 \ 
where = In 2/9, a 2 = In 2/25, and 5 is an amplitude factor. By choosing a small 
5 = 0.001, the Euler equations are practically linearized. During the computation, 
the NRBC of Type-I described above is enforced at all boundaries. Fig. 2 illustrates 
the density contours at different time steps. It also shows that the simple NRBC is 
effective. Comparison of perturbed density distributions along the centerline y = 0 
with the exact solutions at t = 60 and t = 100 are demonstrated in Plots (a) and 
(b) respectively. The numerical results agree well with the theoretical ones, with 
no visible dispersion error. 

3.2. Linear Instability of a Free Shear Layer and Vortex Roll-Up. In 

the 2nd example [5], the inviscid linear and nonlinear instability properties of a 
free shear layer is studied. The numerical results are compared with linear results 
obtained by the normal mode approach. 

The background mean flow consists of a fast (supersonic) stream on the top 
half domain and a slow (subsonic) stream at the bottom half, with the two parallel 
streams connected by a shear layer of hyperbolic tangent shape. The nondimen- 
sional flow states are given as: 

U\ = 1, Vi = 0, pi = 1/3.15, pi = l, 

U 2 = .7391304, V 2 = 0, p 2 = 1/3.15, p 2 = 0.5405405. 

The computational domain spans between 0 < x < 300 and — 10 < y < 10, with 
a uniform structured grid of 601 x 101 (1201 x 201 for fine grid), with time step 
size At = 0.15. The computation is carried out until t = 390 when the spatial 
instability is fully developed. A small harmonic perturbation (amplitude 5 = 0.001) 
at the most unstable frequency / = 0.062 is enforced at the inlet boundary. Fig. 3 
shows the power spectrum P n of the computed u ' (perturbation of u) in natural 
log scale at the streamwise stations x = 50, 100, 150, and 250, which correspond 
to about 3.5, 7 wavelengths and so on. At x = 50, there is a clearly discernible 
peak centered at about the forcing frequency / = 0.062. At about x = 100, 
second and third harmonics start to emerge. Further downstream, more harmonics 
appear and eventually the fundamental saturates. The streamwise evolution of 
the disturbance amplitude along the horizontal centerline y — 0 is illustrated in 
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Figure 2. Propagation of an acoustic pulse with a vortic- 
ity /entropy disturbance, and comparison to exact solution. 


Fig 4. Numerical results for both the coarse and fine grids are presented and 
compared to the theoretical linear growth. It is seen that for 0 < x < 100 (about 7 
wavelengths) the coarse-grid and fine-grid computations yield good agreement both 
between themselves and with the linear result. Further downstream, nonlinear 
effects become important, the growth rate is reduced from the linear value, and 
ultimately the fundamental saturates. Both numerical results agree reasonably 
well with each other also in this nonlinear region. Fig. 5 shows the streamwise 
evolution of the disturbance phase. The numerical results for both the coarse 
and finer grids are compared to the corresponding result from linear theory. The 
agreement is surprisingly good until well into the nonlinear region. Fig. 6 compares 
the normalized |^'| profile across the shear layer flow with the eigenfunction from 
linear stability theory at the streamwise station x = 100, which is located towards 
the end of the linear growth region. The agreement is excellent. The phase variation 
across the shear layer of the disturbance at the same station is depicted in Fig. 7. 
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Figure 3. Power spectra at x = 050, 100, 150, and 250 with 
forcing at the most unstable frequency according to linear theory; 
(coarse grid) 



Figure 4. Streamwise evolution of disturbance amplitude with 
forcing at the most unstable frequency according to linear theory 
(/ = .062). Squares: total u rmsy coarse grid; triangles: total u rmsy 
finer grid; circles: u rms at forcing frequency, coarse grid; crosses: 
u rms at forcing frequency, finer grid; solid line: linear growth. 


Fig. 8 shows the contours of the flow variables for the finer-grid computation. This 
figure clearly demonstrates the effectiveness of the NRBC’s at the top, bottom and 
outlet. It should be emphasized that the domain shown in the figures is exactly 
the computational domain, no buffer zones, cut-offs or other numerical fixes are 
applied. 


3.3. Mach radiation from a supersonic axisymmetric jet. Another in- 
teresting linear or quasi-linear wave phenomenon is the Mach radiation from a fully 
expanded supersonic jet. In the following examples, a perturbation (or forcing) is 
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Figure 5. Streamwise evolution of disturbance phase with forcing 
at the most unstable frequency according to linear theory (f=.062). 
Circles: coarse grid; crosses: finer grid; solid line: linear result. 



Figure 6 . Transverse mode shape at x = 100 with forcing at most 
unstable frequency according to linear theory, finer grid. Squares: 
total u rms ; circles: u rms at forcing frequency; solid line: linear 
eigenfunction (modulus). 


provided by a right-hand-side source term in the energy equation, the 4th compo- 
nent of Eq. (1), located at the origin (0, 0) inside the jet core: 

S 

exp [~B(x 2 + y 2 )\ cos (cot), 

7-1 

where 5 = 0.001 is a small number, c o = 27 iSt (St being the Strouhal number), 
and the constant B = 8. In this test example [6], a fully expanded supersonic 
jet with Mach number Mj = 2.0 is considered. The computational domain spans 
between 0 < x < 33 D and 0 < y < 19 D, with a non-uniform structured grid 
of 350 x 281 nodes, where D is the diameter of the jet at the nozzle exit. More 
grid points are packed around the shear layer. The computation of the unsteady 
jet flow is carried out to t = 100 when the spatial instability is fully developed. 
Fig. 9 illustrates the isobars and v- velocity contours in the near and intermediate 
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Figure 7. Transverse phase variation at x = 100 with forcing 
at most unstable frequency / = .062 according to linear theory, 
(fine grid). Circles: u rms at forcing frequency; solid line: linear 
eigenfunction (phase). 


field. The Mach radiation flow pattern agrees qualitatively with experimental [7] 
and other computational [8] results. 

4. 2-D Nonlinear aeroacoustics Test problems 

In this section, several examples of non-linear aeroacoustics involving shocks 
are demonstrated to show the capability of the CE/SE scheme. 

4.1. Multiple Interaction of a Strong Vortex and Shocks. Vortex-shock 
interaction has been considered a ’difficult’ problem for CAA , since most of the 
CAA schemes handle weak acoustic waves and srong socks simultaneously. How- 
ever, the CE/SE Euler scheme can solve the problem in an effortless way, without 
any numerical trick or fix [5]. As seen in Fig 10, a uniform structured grid of 
401 x 101 nodes is employed in the rectangular domain with Ax = Ay = 1. The 
inflow boundary condition is given as a supersonic flow with a Mach number of 2.9: 

= 2.9, v 0 = 0, po = 1/1.4, po = 1. 

The boundary condition at the top is an inclined flow: 

u t = 2.6193, v t = -0.50632, p t = 1.5282, p t = 1.7000. 

The outflow boundary condition is the type-II NRBC and the bottom is a solid 
reflecting wall. Then, a steady oblique shock is formed with 29° inclination and 
reflected at the bottom wall. The flow with shocks is pre-calculated until a steady 
state is reached. It is then used as the background mean flow for further computa- 
tion. 

At t = 0, a strong Lamb’s vortex is placed at (22,60) (Fig. 10). With At = 0.2, 
900 time steps were run. Fig. 10 demonstrates the shock- vortex interaction at the 
different times t = 2, 20, 38, 56, 74, 92, 110, 128, 146, and 180, and shows how the 
nonlinear acoustic waves are generated, and how they pass through the shocks and 
convect downstream. 
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Figure 8. Flow variable contours at t — 390, fine grid 


4.2. Screech tone noise of an axisymmetric underexpanded jet. Screech 
tone noise is a typical nonlinear aeroacoustic phenomenon. In nature, the screech 
noise is caused by an acoustic feedback loop between the jet nozzle exit lip (or 
nearby solid surfaces) and the shock cell structure. At the low supersonic jet Mach 
number Mj = 1.19, the overall motion in the experiment [8] is in an axisymmetric 
mode, and the use of an axisymmetric 2-D scheme is appropriate. 

The Navier-Stokes CE/SE solver with LES (large eddy simulation) is applied to 
a triangulated grid of 88300 cells [10]. The geometry of the computational domain 
is shown in Fig 11. Both Type I and II NRBC’s are applied appropriately to 
all boundaries except at the nozzle solid walls where no-slip solid wall boundary 
condition is employed, and at the nozzle exit where the following inflow condition 
is applied: 

i 

7 -l 


Pe =P0 


'l + |(7-l )M] 
1 + |(7 — 1) 


p e — 0.5(7 + 1)P0 


1 + !( 7 _1) M 2- 

1 + 1(7— 1) 
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Figure 9. Mach radiation of a supersonic jet at Mj = 2.0 and 
Strouhal number St = 0.2, domain size 33Dxl9D, non-uniform 
grid. 



“ e -(( 7 + i)) 1/2 ’ We -°* 

Initially the entire field is set to the quiescent ambient flow: 


Po = 1, uo — 0, v 0 = 0, po = 

7 


In the numerical simulation, the initial impact of the boundary condition at the 
nozzle exit will stimulate the jet shear layer and trigger the feedback loop to form 
the self-sustaining oscillation that generates screech waves. Fig. 13 illustrates the 
radiating screech waves at the time level of 410,000 steps. Since there is no forcing 
at all, these waves are a sure sign of the presence of a sustainable self-excited oscil- 
lation. In processing the numerical data for this figure, the very high level isobar 
contours, corresponding to hydrodynamic waves around the jet core area are ‘cut 
off’ and the ‘colors’ are appropriately adjusted so that the acoustic waves are clearly 
displayed. The screech wavelength is about 1.6D, well in line with experimental 
results (e.g. [8]). The shock cell structure in Fig. 13 is comparable to the exper- 
imental Schlierens in Fig. 12. Fig. 14 shows the PSD (power spectrum density) 
at the location (x = 2.0, y = 6.0). It is seen that the spikes at the fundamental 
frequency and the subharmonics shoot out against the background noise, despite 
the low resolution. The fundamental frequency is estimated to be about 8, 661 Hz 
in this case, comparable to the experimental case - 8,525 Hz. 
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isobars for a vortex passing through shocks, with 
acoustic waves generated. 


Figure 10. Multiple interactions of a strong vortex and shocks 
producing nonlinear acoustic waves. 


4.3. Transonic resonance of a convergent- diver gent (C-D) nozzle. The 

transonic resonance phenomenon of a C-D nozzle is a newly discovered jet noise 
[11] and is believed to be a precursor of such jet noises as the screech tones. Fig 15 
is a sketch of the C-D nozzle transonic resonance problem. When the flow is at 
transonic speed, a shock Mach disc is formed at the nozzle throat, the disc vibrates 
like a diaphragm as a result of the self-sustaining oscillation in the flow, and triggers 
the quarter wave-length resonance of the divergent part of the nozzle, as if it were a 
circular duct with one end close. Fig. 16 shows that the computed frequencies agree 
very well with the experimental ones, including a stage jump where the resonance 
switches from fundamental frequency to the third harmonic. Furthermore, the two 
plots in Fig. 17 demonstrate the r.m.s. pressure fluctuation along the centerline 
of the divergent portion of the nozzle, showing clearly the fundamental frequency 
resonance (bottom plot) and third harmonic resonance (top plot). The numerical 
work thus helps to confirm the theory on the transonic resonance. The numerical 
work is detailed in [11]. 

4.4. the 3-D screech noise from a circular jet. The present problem is 
a 3-D version of the jet screech problem in §4.2. In order that the flow and and 
noise are truly three dimensional, the jet Mach number is set at Mj = 1.42 [9]. 
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Figure 11. Geometry of the computational domain. 



Figure 12. Schlieren picture from experiment (Panda), showing 
shock cells. 

Approximately 3.7 million hexahedral cells are employed in the numerical solution. 
For a huge grid like this, parallel computation becomes necessary. A Linux PC 
cluster of 60 processors is used. Fig. 18 shows that the computed shock cell structure 
agrees well with the experimental one. For example, the shock cell spacing is about 
1.28.D, D being the jet diameter at the nozzle exit. The radiated screech waves are 
in the flapping B mode, as shown in Fig. 19. 

5. The Airframe Noise Problems 

Airframe noise is an important noise source in aeroacoustics. Several airframe 
tone noise problems are considered in this section to demonstrate the capability of 
the CE/SE scheme. 

5.1. The subsonic cavity noise problem. In this problem[12], a M — 0.8 
subsonic flow passes over a cavity of aspect ratio 6.5. Tonal oscillations occur in 
a feedback cycle in which the vortices shed from the upstream lip of the cavity 
convect downstream and impinge on the other lip edge, generating acoustic waves 
that in turn propagate upstream to excite new vortices (Fig 20). 

In the computational domain, 42,000 triangular elements are used in the un- 
structured grid, which is made from structured rectangular cells. Figure 21 depicts 
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isobars at time step 410,000 



Figure 13. Screech wave radiation and shock cell structure (nu- 
merical Schlierens). 
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Figure 14. PSD at (2,6); x-axis: Strouhal number St, St = 1 is 
equivalent to 13504 Hz. 


with isobars where the acoustic waves are generated and propagate in a series of 
snapshots (1-12) in the near field of the cavity. Each snapshot is 3.6 (720 steps) unit 
apart in time. Initially, the flow conditions are set to be identical to the quiescent 
ambient flow. The boundary conditions are the appropriate types of NRBC except 
at the cavity walls, where no slip wall boundary condition is applied. No visible 
reflections are observed at the non-reflecting boundaries. At the inflow boundary, 
upstream propagating waves are well absorbed within 3-5 cells, without contam- 
inating the interior domain. From an animation of the solution, the near field 
acoustic wave structure appears to be complicated and chaotic. 
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Figure 15. Transonic resonance of a C-D nozzle. 



Figure 16. Frequency vs. jet Mach number Mj, showing fre- 
quency jumps (staging). 


5.2. Flat plate trailing edge noise. In this last example, the noise gen- 
erated by the von Karman vortex street downstream of a rectangular blunt slab 
is computed [13]. The flow Mach number is M — 0.3. Figure 22 illustrates the 
details of the computational domain. Initially, the entire field is set at the ambient 
flow. Then the inflow of M — 0.3 is imposed at the left boundaries. At the solid 
slab walls, no-slip condition is applied. The computation is totally carried out for 
390, 000 time steps. Time history at the data point as shown in Fig. 22 is recorded 
for FFT analysis. The power spectrum density (PSD) is shown in Fig. 23. As be- 
fore, the x— axis denotes the reduced frequency - the Strouhal number5t. St = 1 
is equivalent to a frequency of 13, 504 Hz if the thickness of the slab is assumed to 
be one inch. The frequency at the lowest tone spike is 791 Hz, almost identical to 
Heinemann’s experimental result [14]. 


NAS A/CR— 2003 -2 12388 


14 



. <b) 

9 

0 9 

f \ 



i 

\ 0 \ 

© 

1 i \ 

k 

-Vj 

1 1 1 1 1 


0 " 

<z> 

O 

i i i i 

i i i i i 



Figure 17. r.m.s. pressure fluctuation along the nozzle centerline. 




Figure 18. Comparison of experimental and numerical schlierens, 
showing shock cell structure. Note that the experimental result is 
a time- averaged one and the numerical one is instantaneous. 

6. Concluding remarks 

Through numerous numerical examples, the capability of the recent CE/SE 
schemes for aeroacoustics computations is demonstrated. The scheme features: 

(1) naturally adapted to unstructured grid, 

(2) high resolution, low dispersion, and low dissipation, despite its nominal 2nd 
order accuracy (in space and time), 

(3) robust, treats the ’difficult’ problems in a simple, effortless way, particularly 
appropriate for near field, non-linear aeroacoustics, 

(4) the novel NRBC is simple, effective and truly multi-dimensional. 

The CE/SE scheme thus is a viable tool for CFD and CAA, and in particular, 
for nonlinear aeroacoustics. 
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